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1. Introduction 

Simulations in lattice QCD have advanced remarkably in the past couple of years reaching the 
physical pion mass. The theoretical and algorithmic improvements, combined with the tremendous 
increase in computational power, have made ab initio calculations of key observables on hadron 
structure in the chiral regime feasible enabling comparison with experiment. Form factors and 
generalized parton distribution functions (GPDs) can be obtained from the generalized form factors 
in certain limiting cases. GPDs provide detailed information on the internal structure of hadrons 
in terms of both the longitudinal momentum fraction and the total momentum transfer squared. 
Beyond the information that the form factors yield, such as size, magnetization and shape, GPDs 
encode additional information, relevant for experimental investigations, such as the decomposition 
of the total hadron spin into angular momentum and spin carried by quarks and gluons. GPDs 
are single particle matrix elements of the light-cone operator [jlj Qj, which can be expanded in 
terms of local twist-two operators 0p { ^ 2 ' = xjff iD^-iD ^ V J - Lattice QCD allows 
us to extract hadron matrix elements for the twist-2 operators, which can be expressed in terms of 
generalized form factors. 

In order to compare hadron matrix elements of these local operators to experiment one needs 
to renormalize them. The aim of this paper is to calculate non-perturbatively the renormalization 
factors of the above twist-two fermion operators within the twisted mass formulation. We show 
that, although the lattice spacings considered in this work are smaller than 1 fm, 0(a ) terms are 
non-negligible and introduce significantly larger uncairtainties than statistical errors. We therefore 
compute the &(a 2 ) terms perturbatively and subtract them from the non-perturbative results. This 
subtraction suppresses lattice artifacts considerably depending on the operator under study and 
leads to a more accurate determination of the renormalization constants [J3J, |J. 



2. Formulation 

For the gauge fields we use the tree-level Symanzik improved gauge action [|5]|, which includes 
besides the plaquette term also rectangular (1x2) Wilson loops. The fermionic action for two 
degenerate flavors of quarks in twisted mass QCD is given by 

Sf = a A YjC{x){Dw[U] +m + im5^)x(x) (2.1) 

X 

with T 3 the Pauli matrix, /Iq the bare twisted mass and D w the massless Wilson-Dirac operator. 
Maximally twisted Wilson quarks are obtained by setting the untwisted bare quark mass thq to 
its critical value m cr , while the twisted quark mass parameter is kept non-vanishing in order 



to give the light quarks their mass. In Eq. (|2.l|) the quark fields % are in the so-called "twisted 



basis". The "physical basis" is obtained for maximal twist by the simple transformations \jf( x ) = 
exp (f y 5 T 3 ) x(x), xj/(x)=x(x)exp (f y 5 T 3 ). 

Here we consider only the vector and axial twist-two operators up to one-derivative, Zy, Za, 
Zdv> Zda (symmetrized over two Lorentz indices and traceless), which are given in the twisted 
basis as follows: 
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Vrsr^V = 1 _ iwrsr^D v} t 2 \i/ =i 

^ft* 8 * = { -vr5Yn*W o = 2 ^Dv v}n =lr(fi^v}^z = < -WsriM^vj^V = 2 

y/fyrV a = 3 [vr{M D v}^V « = 3 

If** o=l fvfy^v^V o = l 

-TOtV o = 2 ^ v> =Z757{m < Ov}T' , Z = < -fr^D v} T l \lf a = 2 (2.2) 

WsYn^V o = 3 [w5 7{m^v}tV o = 3 

In a massless renormalization scheme the renormalization constants are defined in the chiral limit, 
where isospin symmetry is exact. Hence, the same value for Z is obtained independently of the 
value of the isospin index a and therefore we drop the a index from here on. However, one must 
note that, for instance, the physical YY{n D v ]l y is renormalized with Z^a, while \j/Y{n D V }? 3 Y 
requires the Zdv, which differ from each other even in the chiral limit. The one-derivative op- 
erators fall into different irreducible representations of the hypercubic group, depending on the 
choice of indices. Hence, we distinguish between ^bvi (^dai) = ^bv (<^da) with /I = v and 

^DV2 (^DA2) = ^DV (^Da) with /i ^ V. 

2.1 Renormalization Condition 

The renormalization constants are computed both perturbatively and non-perturbatively in 
the Rf-MOM scheme at various renormalization scales. We translate them to the MS-scheme 
at (2 GeV) 2 using a conversion factor computed in perturbation theory to €?(g 6 ) as described in 
Section 3. The Z-factors are determined by imposing the following conditions: 



Z„ = — Tr 
' 12 



- iL ^{s\ P ))-^ 



P- 



P 2 =H 2 



= 1, (23) 



where pL is the renormalization scale, Si and Tl correspond to the perturbative or non-perturbative 
results and p v \ is the tree-level expression of the operator under study. The trace is taken over 
spin and color indices, and the conditions are imposed in the massless theory. 

2.2 Perturbative procedure 

Our calculation for the Z-factors is performed in 1-loop perturbation theory to iff (a 2 ). The 
order a 2 -terms can be subtracted from non-perturbative estimates, and they can eliminate possible 
large lattice artifacts. There are many difficulties when extracting powers of the lattice spacing 
from our expressions, since there appear singularities encountered at tff(a 2 ), that persist even up to 
6 dimensions (integral convergence in 7-d), making their extraction more delicate. In addition to 
that, there appear Lorentz non-invariant contributions in <^(a 2 ) -terms, such as L^/^/p 2 , where p 
is the external momentum; as a consequence, the Z-factors also depend on such terms. 

For all our perturbative results we employ a Wilson-type fermion action (Wilson/clover/twisted 
mass), with non-zero bare mass, m. For the renormalization of the fermion field and the local bi- 
linears we also have a finite twisted mass parameter, jUo, so we can explore the mass dependence. 
For gluons we use Symanzik improved actions (Plaquette, Tree-level Symanzik, Iwasaki, TILW, 
DBW2) [g]. The expressions for the matrix elements and the Z-factors are given in a general co- 
variant gauge, and their dependence on the coupling constant, the external momentum, the masses 
and the clover parameter csw is shown explicitly. The Feynman diagrams involved in the compu- 
tation of the various Z-factors are illustrated in Fig. [l[ Here we do not show any expressions for the 
matrix elements of the Green's functions, since they are far too lengthy. As an example we show 
the G(a 2 ) terms that can improve the non-perturbative estimate of Z q once they are subtracted. For 
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Figure 1: One-loop diagrams contributing to the correction of the amputated Green's functions of the propa- 
gator (la, lb), local bilinears (2) and one-derivative operators (3a-3d). A wavy (solid) line represents gluons 
(fermions). A cross denotes an insertion of the operator under study. 



the special choices: csw = 0, r = 1 (Wilson parameter), A = (Landau gauge), nig = 0, plo = 0, 
and for tree-level Symanzik gluons, Z q can be corrected to 6{cr) as follows: 



(2.4) 



Its most general expression is far too lengthy to be included in paper form; it is provided, along 
with the rest of our results for the Z-factors, in electronic form in Ref. [Q]. 

2.3 Non-perturbative calculation 

For each operator we define a bare vertex function given by 

G(P) = °y £ e-^-yHuix^jrMdimy)) , (2.5) 

x ,y,z,z' 

where p is a momentum allowed by the boundary conditions, V is the lattice volume, and the 
gauge average is performed over gauge-fixed configurations. The form of J* {z,z') depends on 
the operator under study, for example J? {z,z l )=8 Z)Z iy il would correspond to the local vector cur- 
rent. In the literature there are two main approaches that have been employed for the evaluation of 



Eq. (2.5). The first approach relies on translation invariance to shift the coordinates of the correla- 



tors in Eq. (2.5) to position z=0 [|7|]. Having shifted to z=0 allows one to calculate the amputated 



vertex function for a given operator Jf for any momentum with one inversion per quark flavor. 
In this work we explore the second approach, introduced in Ref. which uses directly Eq. ( ^3| ) 
without employing translation invariance. One must now use a source that is momentum dependent 
but can couple to any operator. For twisted mass fermions, with twelve inversions one can extract 
the vertex function for a single momentum. The advantage of this approach is a high statistical 
accuracy and the evaluation of the vertex for any operator including extended operators at no sig- 
nificant additional computational cost. We fix to Landau gauge using a stochastic over-relaxation 
algorithm [pj]. 



3. Results 



We perform the non-perturbative calculation of renormalization constants for three values of 
the lattice spacing, t?=0.089 fm, 0.070 fm, 0.056 fm, corresponding to j3 = 3.9, 4.05 and 4.20 
respectively. In Tables I and II of Ref. [Q] we summarize the various parameters that we used 
in our simulations. We have tested finite volume effects and pion mass dependence; both effects 
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are within the small statistical errors for the operators considered here. Chiral extrapolations are 
necessary to obtain the renormalization factors in the chiral limit. Since the dependence on the pion 
mass is insignificant, even if we allow a slope and perform a linear extrapolation to our data, this 
is consistent with zero; therefore the renormalization constants are computed at one quark mass. 
Figures ||||3] demonstrate the effect of subtraction at two j8 values for the local and one-derivative 
vector/axial Z-factors, as a function of the renormalization scale (in lattice units). Zy and Za 
are scale independent, thus we obtain a very good plateau upon subtraction of 0(a 2 ) effects. To 
identify a plateau for Zdv and Zda we need to convert to MS and evolve to a reference scale. 
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Figure 2: Renormalization scale dependence for Zy, Za at j5 — 3.9, m K = 0.430 GeV (left panel) and 
/3 = 4.20, m K = 0.476 GeV (right panel) (Open points: unsubtracted, filled points: subtracted). 




Figure 3: Renormalization scale dependence for Zdv, Zda (RI'-MOM scheme before evolving at a reference 
scale) at j3 = 3.9, m x = 0.430 GeV (left panel) and j3 = 4.20, m n = 0.476 GeV (right panel) (Open points: 
unsubtracted, filled points: subtracted). 

• Conversion to MS: The passage to the continuum MS -scheme is accomplished through use of 
a conversion factor, which is computed up to 3 loops in perturbation theory. By definition, this 
conversion factor is the same for the one-derivative vector and axial renormalization constant, but 
will differfor the cases Z DV i (Z DA i) and Z DV2 (Z DA 2), that is C D vi = Cdai = z dv/ z dvi > C DV2 = 
Cda2 = z dv/ z dv2- This requirement for different conversion factors results from the fact that the 
Z-factors in the continuum MS-scheme do not depend on the external indices, jj., v (see Eq. (2.5) 



of Ref. [1C]), while the results in the RI'-MOM scheme do depend on /I and v. We also need 
another factor R(2GeV,jJ.) that will bring all Z-factors down to jj. = 2 GeV, for example 



Zf^ (2GeV) = R DY (2GeV,ii) ■ C DV i (jl) ■ Zg wl Qi) 



RI' 



(3.1) 
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Figure 4: Renormalization factors at j3 = 3.9, m K = 0.430 GeV (left panel) and j3 = 4.20, m K = 0.476 GeV 
(right panel) in the MS-scheme at 2 GeV. The lines show extrapolations to a 2 p 2 = within the range p 2 ~ 
15 — 32 (GeV) 2 . (Open points: unsubtracted, filled points: subtracted) 

A "renormalization window" should exist for Aq CD « \i 2 << I /a 2 where perturbation the- 
ory holds and finite-a artifacts are small, leading to scale-independent results (plateau). In practice 
such a condition is hard to satisfy: The upper range of the inequality is extended to (2 — 5) /a 2 lead- 
ing to lattice artifacts in our results that are of &(a 2 p 2 ). Fortunately our perturbative calculations 
allow us to subtract the leading perturbative 0{a 2 ) lattice artifacts which alleviates the problem. To 
remove the remaining 0(a 2 p 2 ) artifacts we extrapolate linearly to a 2 p 2 = as demonstrated in Fig. 
|J The statistical errors are negligible and therefore an estimate of the systematic errors is impor- 
tant. We note that, in general, the evaluation of systematic errors is difficult. The largest systematic 
error comes from the choice of the momentum range to use for the extrapolation to a 2 p 2 = 0. One 
way to estimate this systematic error is to vary the momentum range where we perform the fit. An- 
other approach is to fix a range and then eliminate a given momentum in the fit range and refit. The 
spread of the results about the mean gives an estimate of the systematic error. In the final results we 
give as systematic error the largest one from using these two procedures which is the one obtained 
by modifying the fit range. In order to treat all beta values equally, we fix the momentum range in 
physical units and we thus fit all renormalization constants in the same physical momentum range, 
p 2 ~ 15 — 32 (GeV) 2 . The momentum interval in physical units has bean chosen such as a good 
plateau exists at each j8, as can be seen in Fig. |]. The Gia 1 ) perturbative terms which we subtract, 
decrease as j8 increases, as expected. The momentum range in lattice units at each /3 is rescaled as 
follows: p = 3.9 : a 2 p 2 ~ 3 - 5, J8 = 4.05 : a 2 p 2 ~ 1.9 -3, j3 = 4.20 : a 2 p 2 ~ 1.2 - 2.5. Our re- 
sults for the €?(a 2 ) corrected Z-factors in the MS-scheme at 2 GeV are given in Table [T[ which have 
been obtained by extrapolating linearly in a p . For Zpy and Zda we used the fixed momentum 
range p 2 ~ 15 — 32 (GeV) 2 [Q], while for Zy and Za we used all the data points available, since the 
plateau is good for all momenta. The final results for Zy and Za for a more extended momentum 
range will appear in [Jj]. 



4. Conclusions 

The values of the renormalization factors for the one-derivative twist-2 operators are calculated 



6 



Renormalization constants for one-derivative fermion operators in twisted mass QCD M. Constantinou 



p 


Zv 


Z A 


Zdvi 


ZdV2 


Zdai 


ZdA2 


3.90 
4.05 
4.20 


0.6343(6)(3) 

0.6628(7)(14) 

0.6854(5X13) 


0.7561(6)(5) 
0.7722(6X3) 
0.7870(5X9) 


0.970(34)(26) 
1.033(11)(14) 
1.097(4)(6) 


1.061(23)(29) 
1.131(23)(18) 
1.122(7)(10) 


1.126(22)(78) 

1.157(9)(7) 

1.158(7)(7) 


1 .076(5)(1) 

1.136(5) 

1.165(5)(10) 



Table 1: Renormalization constants in the MS scheme, after extrapolating linearly in a 2 p 2 . The error in the 
first parenthesis is statistical and the one in the second parenthesis is systematic. 



non-perturbatively. The method of choice is to use a momentum dependent source and extract the 
renormalization constants for all the relevant operators, which leads to a very accurate evaluation 
of these renormalization factors using a small ensemble of gauge configurations. We studied the 
quark mass dependence and found that an extrapolation to zero quark mass changes the result by 
about 1 per mille for all the operators we presented here. This is in most cases by an order of 
magnitude smaller than the systematical errors due to lattice artifacts, therefore a calculation at a 
single quark mass suffices. For all the renormalization constants shown here we do not find any 
light quark mass dependence within our small statistical errors. Therefore it suffices to calculate 
renormalization constants at a given quark mass. Despite using lattice spacing smaller than 1 fm, 
0(a 2 ) effects are sizable, thus, we perform a perturbative subtraction of ff(a 2 ) terms. This leads 
to a smoother dependence of the renormalization constants on the momentum values at which they 
are extracted. Residual &{a rp > ) effects are removed by extrapolating to zero. In this way we can 
accurately determine the renormalization constants in the Rf-MOM scheme. In order to compare 
with experiment we convert our values to the MS scheme at a scale of 2 GeV. The systematic errors 
are estimated by ochanging the window of values of the momentum used to extrapolate to a 2 p 2 = 0. 

References 

[1] Xiang-Dong Ji, J. Phys. G24(1998) 1181, [hep-ph/9807358]. 

[2] LHPC Collaboration: Ph. Hagler, J.Negele, D. Renner, W. Schroers, Th. Lippert, K. Schilling, Phys. 
Rev. D68, (2003) 034505, [hep- 1 at / 3 4 1 8] . 

[3] C. Alexandrou, M. Constantinou, T. Korzec, H. Panagopoulos, F. Stylianou, accepted in Phys. Rev. D, 

[arXiv: 1006 . 1920]. 

[4] C. Alexandrou, M. Constantinou, T. Korzec, H. Panagopoulos, F. Stylianou, in preparation. 
[5] P. Weisz, Nucl. Phys. B212 (1983) 1. 

[6] M. Constantinou, V. Lubicz, H. Panagopoulos, F. Stylianou, /HEP 10 (2009) 064, 
[arXiv: 0907 . 0381]. 

[7] ETM Collaboration: M. Constantinou, P. Dimopoulos, R. Frezzotti, G. Herdoiza, K. Jansen, V. 
Lubicz, H. Panagopoulos, G.C. Rossi, S. Simula, F. Stylianou, A. Vladikas, /HEP 08 (2010) 068, 

[arXiv: 1004 .1115]. 

[8] M. Gockeler, R. Horsley, H. Oelrich, H. Perlt, D. Petters, P.E.L. Rakow, A. Schafer, G. Schierholz, A. 
Schiller, Nucl. Phys. B544 (1999) 699, [hep-lat/9807 04 4]. 

[9] Ph. de Forcrand, Nucl. Phys. Proc. Suppl., 9 (1989) 516. 

[10] J. A. Gracey, Nucl. Phys. B667 (2003) 242, [hep-ph/030 61 63]. 



7 



